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ABSTRACT 

We present HST/STIS spectra of the She spiral galaxy NGC 4041 which were used to map 
the velocity field of the gas in its nuclear region. We detect the presence of a compact (r ~ 0"4 ~ 
40 pc), high surface brightness, rotating nuclear disk co-spatial with a nuclear star cluster. The 
disk is characterized by a rotation curve with a peak to peak amplitude of ~ 40 km s“^ and is 
systematically blueshifted by ~ 10—20 km s“^ with respect to the galaxy systemic velocity. With 
the standard assumption of constant mass-to-light ratio and with the nuclear disk inclination 
taken from the outer disk, we find that a dark point mass of (I^q®) x 10^ Mq is needed to 
reproduce the observed rotation curve. However the observed blueshift suggests the possibility 
that the nuclear disk could be dynamically decoupled. Following this line of reasoning we relax 
the standard assumptions and find that the kinematical data can be accounted for by the stellar 
mass provided that either the central mass-to-light ratio is increased by a factor of ~ 2 or that the 
inclination is allowed to vary. This model results in a 3cr upper limit of 6 x 10® Mq on the mass of 
any nuclear black hole. Overall, our analysis only allows us to set an upper limit of 2 x 10^ Mq on 
the mass of the nuclear BH. If this upper limit is taken in conjunction with an estimated bulge 
B magnitude of —17.7 and with a central stellar velocity dispersion of ~ 95 km s“^, then these 
results are not inconsistent with both the MeH-T^sph and the Mbh-c* correlations. Constraints 
on BH masses in spiral galaxies of types as late as Sbe are still very scarce and therefore the 
present result adds an important new datapoint to our understanding of BH demography. 

Subject headings: black hole physics — galaxies: individual (NGC 4041) — galaxies: kinematics and 
dynamics — galaxies: nuclei — galaxies: spiral 
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1. Introduction 

It has long been suspected that the most lu¬ 
minous AGN are powered by accretion of mat- 
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ter onto massive black holes (hereafter BH; e.g. 
Lynden-Bell 1969). This belief, combined with the 
observed evolution of the space-density of AGN 
(Soltan 1982; Chokshi & Turner 1992; Marconi & 
Salvati 2001) and the high incidence of low lu¬ 
minosity nuclear activity in nearby galaxies (Ho, 
Filippenko, & Sargent 1997), implies that a signif¬ 
icant fraction of luminous galaxies must host black 
holes of mass 10® - 10^° M©. 

It is now clear that a large fraction of hot 
spheroids (E-SO) contains a BH (Harms et al. 1994; 
Kormendy & Richstone 1995; Kormendy et al. 
1996; Macchetto et al. 1997; van der Marel & van 
den Bosch 1998; Bower et al. 1998; Marconi et al. 
2001) with mass proportional to the mass (or lumi¬ 
nosity) of the host spheroid (MBn/AIsph ~ 0.001 
e.g. Merritt & Ferrarese 2001). Recently Ferrarese 
& Merritt (2000) and Gebhardt et al. (2000) have 
shown that a tighter correlation holds between the 
BH mass and the velocity dispersion of the bulge. 
Clearly, any correlation of black hole and spheroid 
properties would have important implications for 
theories of galaxy formation in general, and bulge 
formation in particular. However, to date, there 
are very few secure BH measurements or upper 
limits in spiral galaxies even though we know that 
AGNs are common in such systems (Maiolino & 
Rieke 1995). In total, there are 37 secure BH 
detections according to Kormendy & Gebhardt 
(2001) or just 22 according to Merritt & Ferrarese 
(2002), depending on the dehnition of “secure”. 
Only ~ 20% of these BH detections (7/37 or 4/22, 
respectively) are in galaxy types later than SO, and 
only 3 in Sbc’s and later (the Milky Way, Genzel 
et al. 2000; NGC 4258, Miyoshi et al. 1995; NGC 
4945, Greenhill, Moran, & Herrnstein 1997). It is 
therefore important to directly establish how com¬ 
mon are BHs in spiral galaxies and if they follow 
the same Msu-Msph, Mbh-ct* correlations as El¬ 
liptical galaxies. 

This can be achieved only with a comprehensive 
survey for BHs that covers quiescent and active 
spiral galaxies of all Hubble types. Such a survey 
would pin down the mass function and space den¬ 
sity of BHs, and their connection with host galaxy 
properties (e.g., bulge mass, disk mass etc). 

To detect BHs one requires spectral informa¬ 
tion at the highest possible angular resolution: 
the “sphere of influence” (Bahcall & Wolf 1976) 
of BHs are typically < 1" even in the most nearby 


galaxies. Nuclear absorption-line spectra can be 
used to demonstrate the presence of a BH (Ko¬ 
rmendy & Richstone 1995; Richstone 1998; van 
der Marel et al. 1998b), but the interpretation 
of the data is complex because it involves stellar- 
dynamical models that have many degrees of free¬ 
dom - that can be pinned down only when data 
of very high S/N are available (Binney & Mamon 
1982; Statler 1987; Merritt 1997; Binney & Mer- 
rifield 1998). Radio-frequency measurements of 
masers in disks around BHs provide some of the 
most spectacular evidence for BHs, but have the 
disadvantage that only a small fraction of the disks 
will be inclined such that their maser emission 
is directed toward us (Braatz, Wilson, & Henkel 
1997). Studies of ordinary optical emission lines 
from gas disks provide an alternative and rela¬ 
tively simple method to detect BHs. 

HST studies have discovered many cases of such 
gas disks in early-type galaxies (M87, Ford et al. 
1994; NGG 4261, Jaffe et al. 1996; NGC 5322, 
Carollo et al. 1997; Gen A, Schreier et al. 1998) 
and have demonstrated that both their rotation 
curves and line profiles are consistent with thin 
disks in Keplerian motion (Ferrarese et al. 1996; 
Macchetto et al. 1997; Ford et al. 1998; van der 
Marel & van den Bosch 1998; Bower et al. 1998; 
Marconi et al. 2001). 

In early-type galaxies there are still worrying 
issues about the dynamical conhguration of nu¬ 
clear gas (e.g., misalignment with the major axis, 
irregular structure etc). By contrast nuclear gas 
in relatively quiescent spirals is believed to be or¬ 
ganized into well dehned rotating disks seen in op¬ 
tical line images (e.g., M81, Devereux, Ford, & Ja¬ 
coby 1997). Ho et al. (2002) recently found that 
the majority of spiral galaxies in their survey has 
irregular velocity fields in the nuclear gas, not well 
suited for kinematical analysis. Still, 25% of the 
galaxies where Ha emission was detected all the 
way to the center have velocity curves consistent 
with circular rotation and the galaxies with more 
complicated velocity curves can also be useful for 
BH mass measurement after detailed analysis of 
the spectra. Indeed, even in the most powerful 
Seyfert nuclei such as NGG 4151, where the gas is 
known to be interacting with radio ejecta, it may 
be possible to get the mass of the BH from spa¬ 
tially resolved HST spectroscopy by careful analy¬ 
sis of the velocity Held to separate the underlying 
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quiescently rotating disk gas from that disturbed 
by the jets (Winge et al. 1998). 

Prompted by these considerations, we have un¬ 
dertaken a spectroscopic survey of 54 spirals using 
STIS on the Hubble Space Telescope. Our sam¬ 
ple was extracted from a comprehensive ground- 
based study by Axon et al. who obtained Ha and 
[N H] rotation curves at a seeing-limited resolution 
of 1", of 128 Sb, SBb, Sc, and SBc spiral galax¬ 
ies from RC3. By restricting ourselves to galaxies 
with recession velocities V < 2000 km/s, we ob¬ 
tained a volume-limited sample of 54 spirals that 
are known to have nuclear gas disks and span wide 
ranges in bulge mass and concentration. The sys¬ 
temic velocity cut-off was chosen so that we can 
probe close to the nuclei of these galaxies, and de¬ 
tect even lower-mass black holes. The frequency 
of AGN in our sample is typical of that found in 
other surveys of nearby spirals, with comparable 
numbers of weak nuclear radio sources and LIN¬ 
ERS. The sample is described in detail by Axon 
et al. (paper in preparation). 

This paper presents the observations of NGC 
4041, the first object observed, and a detailed de¬ 
scription of the analysis and modeling techniques 
which will be applied to the other galaxies in the 
sample. From the Lyon/Meudon Extragalactic 
Database (LEDA^^), NGC 4041 is classified as a 
Sbc spiral galaxy with no detected AGN activ¬ 
ity. Its average heliocentric radial velocity from 
radio measurements is 1227 ± 9 km s“^ becoming 
~ 1480 km s“^ after correction for Local Group 
infall onto Virgo. With Ho=75 km s“^ Mpc“^ 
this corresponds to a distance of ~ 19.5 Mpc and 
to a scale of 95 pc/". 

The outline of the paper is as follows. In §2 
we present the adopted observational strategy and 
data reduction techniques. In §3 we present the 
rotation curves of the ionized gas and the broad¬ 
band images of the nuclear region of the galaxy. In 
§4.1 we derive the stellar luminosity density from 
the observed surface brightness distribution which 
is then used in the model fitting of the kinematical 
data. All the details of the inversion procedure are 
described in Appendix A. The model fitting of the 
kinematical data is described in §4.2 and all the 
details of the model computation are described in 
Appendix B. In particular, §4.2.2 describes the 

^^http: //leda.univ-lyonl.fr 
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Fig. 1.— Slit positions overlaid on the acquisition 
image. The 0,0 position is the position of the tar¬ 
get derived from the STIS ACQ procedure. The 
white cross is the kinematic center derived from 
the fitting of the rotation curves (see §4.2). 

standard approach, while in §4.2.3 an alternative 
approach is considered. In §5 we discuss the effects 
of our assumptions on the derived value of the BH 
mass for which, in the present case, only an upper 
limit can be set. We then compare this upper 
limit with the MeH-Tsph and Mbh-o’* correlations. 
Finally, our conclusions are presented in §6. 

2. Observations and Data Reduction 
2.1. STIS observations 

NGC 4041 was observed with STIS on the HST 
in 1999 July 7. An ~ 5" x 5" acquisition image 
was obtained with the F28X50LP filter and the 
galaxy nucleus, present within the field of view, 
was subsequently centered and re-imaged follow¬ 
ing the ACQ procedure. The exposure time of the 
acquisition images was 120 s. 

The observational strategy consisted in obtain¬ 
ing spectra at three parallel positions with the cen¬ 
tral slit centered on the nucleus and the flanking 
ones at a distance of 072. The slit positions are 
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overlaid on the acquisition image in Figure 1 and 
their position angle is 43°. At each slit position 
we obtained two spectra with the G750M grating 
centered at Ha, with the second spectrum shifted 
along the slit by an integer number of detector 
pixels in order to remove cosmic-ray hits and hot 
pixels. The nuclear spectrum (NUC) was obtained 
with the O'.'l slit and no binning of the detector 
pixels, yielding a spatial scale of O'.'0507/pix along 
the slit, a dispersion per pixel of AA = 0.554 A 
and a spectral resolution of 7^ = A/(2AA) ~ 6000. 
The off-nuclear spectra (POSl and POS2) were 
obtained with the 0'.'2 slit and 2x2 on-chip binning 
of the detector pixels, yielding O'.'lOl/pix along the 
slit, 1.108 A/pix along the dispersion direction and 
TZ ~ 3000. Total exposure times were 950 s for the 
NUC position and 420 s and 500 s for POSl and 
POS2, respectively. 

The acquisition images were flat-fielded, re¬ 
aligned and co-added in order to improve the 
signal-to-noise ratio. The pixel scale is 070507. 
The flux calibration was first obtained using the 
PHOTFLAM header in the image. We subse¬ 
quently applied a color correction to convert to 
Johnson R magnitudes. In order to do this, we 
used the average spectra of Sb and Sc spiral galax¬ 
ies from Kinney et al. (1996) as spectral templates. 

The raw spectra were reprocessed through the 
calstis pipeline using the darks obtained daily for 
STIS. Standard pipeline tasks were used to ob¬ 
tain flat-field corrected images. The two exposures 
taken at a given slit position were then realigned 
with a shift along the slit direction (by an integer 
number of pixels) and the pipeline task ocrreject 
was used to reject cosmic rays and hot pixels. Sub¬ 
sequent calibration procedures followed the stan¬ 
dard pipeline reduction described in the STIS In¬ 
strument Handbook (Leitherer et al. 2001), i.e. the 
spectra were wavelength calibrated and corrected 
for 2D distortions. The expected accuracy of the 
wavelength calibration is 0.1 —0.3 pix within a sin¬ 
gle exposure and 0.2 — 0.5 pix among different ex¬ 
posures (Leitherer et al. 2001) which converts into 
~ 3—8 km s“^ (relative) and ~ 5—13 km s“^ (ab¬ 
solute) . The relative error on the wavelength cali¬ 
bration is negligible for the data presented here be¬ 
cause our analysis is restricted to the small detec¬ 
tor region including Ha and [Nil] (AA < lOOA). 

The nominal slit positions obtained as a re¬ 
sult of the STIS ACQ procedure were checked by 


matching the light profiles measured along the slit 
with the synthetic ones derived from the acquisi¬ 
tion image: we collapsed the spectra along the dis¬ 
persion direction and compared the resulting light 
profiles with the ones extracted from the acquisi¬ 
tion image for a given slit position. The agreement 
is good for all slit positions and the center of the 
NUC slit if offset by only ~ 0703 with respect to 
the position of the target determined by the STIS 
ACQ procedure (0,0 position in Figure 1). 

We selected the spectral regions containing the 
lines of interest and subtracted the continuum 
by fitting a linear polynomial row by row along 
the dispersion direction. The continuum sub¬ 
tracted lines were then fitted row by row with 
gaussian functions using the task LONGSLIT in 
the TWODSPEC FIGARO package (Wilkins & 
Axon 1992) and the task specfit in the IRAF^^ 
stsdas package. When the signal-to-noise ratio 
(SNR) was insufficient (faint line) the fitting was 
improved by co-adding two or more pixels along 
the slit direction. 

2.2. WFPC2 images 

WFPC2 images in the F450W (-B), F606W 
(^R) and F814W (~I) filters were retrieved from 
the archive. These images encompass the en¬ 
tire galaxy with the nucleus located in the WF3 
chip. The data were automatically reprocessed 
with the best calibration files available before re¬ 
trieval. Two exposures were performed in each of 
the three filters, in order to remove cosmic rays. 
Warm pixels and cosmic rays were removed (STS¬ 
DAS tasks warmpix and crrej) and mosaicked im¬ 
ages with spatial sampling of 071/pix were ob¬ 
tained using the wmosaic task, which also cor¬ 
rects for the optical distortions in the four WFPC2 
chips. The background was estimated from ar¬ 
eas external to the galaxy where no emission is 
detected. Flux calibration to Vega magnitudes 
was performed using the zero points by Whit¬ 
more (1995). To convert to standard filters in the 
Johnson-Cousins system we estimated the color 
correction using various spectral templates, AOV 
and KOV stars and the Sb and Sc spiral spec- 

^^IRAF is distributed by the National Optical Astronomy 
Observatories, which are operated by the Association of 
Universities for Research in Astronomy, Inc., under coop¬ 
erative agreement with the National Science Foundation. 
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Fig. 2.— From left to right: velocity, FWHM and Surface brightness measured along the slit at POSl, NUC 
and POS2 (from top to bottom). Vertical bars are Icr errors while horizontal bars indicate the size of the 
aperture over which the quantity was measured. In the left panel, the dashed line is the velocity gradient 
measured in the ground based data. The dotted line in the NUC rotation curve panel is the line surface 
brightness along the slit, drawn in order to pinpoint the high surface brightness nuclear disk. The dotted 
lines in the FWHM panels corresponds to instrumental widths. The 0 position along the slit corresponds to 
the position of the nuclear continuum peak. 
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Fig. 3.— Same as previous figure but for the points related to the nuclear disk. The points with the errorbars 
are the quantities derived in the fit which takes into account the presence of the blue wing and in which 
Ha and [Nil] emitting media are constrained to have the same velocity and FWHM. Conversely the points 
without the errorbars are derived with unconstrained, single gaussian fits of Ha (open square) and [N H] 
(crosses). For more details see §3.1.1. 
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Fig. 4.— Fits of the line profiles at the three slit positions. The dotted lines identify the two components 
of the fit. Each component is characterized by the same velocity and FWHM for Ha and [Nil]. The ratio 
between the two [Nil] lines is that fixed by atomic physics. At every slit position the blue component has 
the same velocity and width at each row and their values are determined in the fit of the overall nuclear 
spectrum shown in the upper panels. The numbers in the upper right corners of each panel are the row 
or the range of rows where the fit was performed. The strong narrow lines in the right panels are residual 
cosmic rays which have been excluded from the fit. 


tra from Kinney et al. (1996). The color correc¬ 
tions are the following: / — F814IF ~ —0.1 (Sb, 
KOV) and -0.005 (AOV); i?-F6061F ~ -0.4 (Sb, 
KOV) and -0.05 (Sc,AOV); B - E450IF ~ 0.1 
(Sb, KOV) and 0.0 (Sc,AOV). With these correc¬ 
tions, the differences between colors from Johnson- 
Cousins and HST instrumental magnitudes are 
A(i? - J) ~ -0.3;-0.04, A(B - J) ~ 0.2; 0.0, 
A(i? — i?) ~ 0.5; 0.0 according to the templates 
used. Unless stated otherwise, we have applied 
the color corrections derived from the KOV star 
and Sb spiral templates, which are most suitable 


to the present data. 

2.3. Ground based observations 

NGC 4041 was also observed with NIGS (Baffa 
et al. 2001) at the Telescopio Nazionale Galileo 
(TNG) in 2001 February 12 using the K' filter and 
the small field camera which yields a O'.'IS pixel 
size. The night was not photometric and seeing 
during the observations was ~ O'.'S. Observations 
consisted of several exposures, each with the ob¬ 
ject at a different position on the array. Data re- 
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duction consisted in flat-fielding and sky subtrac¬ 
tion. The frames were then combined into a mo¬ 
saic. In order to flux calibrate the NICS image we 
have used the 2MASS image available on the web. 
The 2MASS K band image is flux calibrated to an 
accuracy of < O.lmag and has a spatial resolution 
of ~ 3'.'5 with a pixel size of 1". We have therefore 
degraded the spatial resolution of the NICS image 
to the 2MASS image and rebinned to pixels. 
We have performed ellipse htting to both images 
and rescaled the NICS image to the flux calibrated 
2MASS data by comparing the light profiles. 

The NICS image and WFPC2 images were then 
realigned by cross correlating common features in 
the nuclear region, because no point sources are 
present in the near infrared image. The main fea¬ 
ture which drives the correlation is the strong nu¬ 
clear peak, however we have checked that, even 
excluding the central peak, the shift between the 
images can be determined with an accuracy of 
±0'.'2 (i.e. ±2 pixels of WFPC2) in both direc¬ 
tions. With this accuracy the nuclear continuum 
peaks present in both optical and infrared images 
are consistent with being at the same position. 
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3. Results 

3.1. Kinematics 

3.1.1. Line fitting procedures 

Line-of-sight velocities, full width at half max¬ 
ima (hereafter FWHM) and surface brightnesses 
along each slit were obtained by htting single gaus- 
sians to Ha and [NII] emission lines in each row 
of the continuum-subtracted 2D spectra. In the 
left panels of Figure 2 we plot the measured ve¬ 
locities. Beyond (/.'b of the nucleus velocities show 
considerable small scale variations likely due to 
local gas motions which do not reflect the mass 
distribution. Therefore we averaged velocities by 
binning the spectra in steps of 1" (10 rows in NUC 
and 5 in POSl,2), while, within 0'.'5 of the nucleus 
velocities are measured along each row to take ad¬ 
vantage of all the spatial information. Similarly, 
the central and right panels in Figure 2 display 
FWHMs and Ha surface brightnesses as measured 
along the slit. 

In Figure 2, the dotted line superimposed on 
the NUC rotation curve represents the Ha sur¬ 
face brightness shown in the right central panel. 


This helps to distinguish the presence of two com¬ 
ponents: an extended one, characterized by a low 
surface brightness (< 10“^^ erg cm“^ s“^ arcsec“^), 
roughly constant along the slit, and nuclear one, 
compact (within —0'.'4 and 0"4), bright and cospa- 
tial with the position of the nuclear continuum 
peak. This might be interpreted as a nuclear disk 
of the same extent as the nuclear stellar cluster 
described in §3.2. As shown in more detail in 
the right hand panel of Figure 3, the emission 
line surface brightness of the nuclear disk is dou¬ 
ble peaked while the nuclear continuum source, 
roughly coincident with the center of rotation, is 
located in-between the two peaks. 

Fitting single unconstrained gaussians is ac¬ 
ceptable for the extended component but it does 
not produce good results for the points in the nu¬ 
clear region. In particular a single gaussian fit 
produces velocities of Ha and [Nil] which differ 
by as much as ~ 30 — 40 km s“^ (see Figure 3). 
This is not a worrisome issue if the amplitude of 
the rotation curve is a few 100 km s“^ but makes 
interpretation of the data uncertain in the present 
case where the amplitude of the nuclear rotation 
curve is only ~ 40 km s“^. A careful analysis of 
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the line profiles in the nuclear region shows that 
they are persistently asymmetric with the presence 
of a blue wing (see Figure 4). 

A fit row by row with two gaussian components, 
the main component and the blue one, with the 
constraint that they have the same velocities and 
widths for Ha and [Nil], shows that, within the 
large uncertainties, the “blue wing” has always the 
same velocity and width. We have then deblended 
the “blue” component in the spectrum obtained 
by co-adding the central 5 rows. The velocity and 
width of the blue component were then used in 
the row by row fit. The constrained fit is good 
and Ha and [N H] now have the same velocity in 
the main component. The measured velocities, 
FWHMs and Ha surface brightnesses in the nu¬ 
clear region are shown in Figure 3 where we also 
plot, as a comparison, the values obtained from 
unconstrained single gaussian fits of Ha and [N H]. 
Given the signal-to-noise of the present data, the 
nature of the blue wing is, as yet, unclear. 

3.1.2. Velocity curves 

The velocity field of the extended component 
shows a quasi-linear gradient of ~ 8 km s~^/” 
which agrees with the velocity gradient (Figure 5) 
measured from ground-based observations at the 
same PA as our STIS spectra (Axon et ah, paper 
in preparation). The velocity gradient measured 
from the ground based rotation curve is shown 
as a straight line in the figures. The large scale 
trend observed from HST matches very well the 
expected gradient from the solid body part of the 
ground based rotation curve although it presents 
structures at small scales. 

Within (y.'5 of the nucleus, the velocity field 
shows a smooth S-shaped curve with a peak to 
peak amplitude of only 40 km s“^ (left panel of 
Figure 3). Overall there is no hint of a steep Ke- 
plerian rise around a point-like mass. 

The central velocities of the nuclear curves are 
systematically offset from the large scale veloc¬ 
ity field by ~ 10 — 20 km s“^ (compare with the 
dashed line which is the solid body part of the 
rotation curve - left panel of Figure 3). The off- 
nuclear slits show essentially the same velocity 
field as the on-nucleus slit. Thus this blueshift 
must be real and not an instrumental artifact gen¬ 
erated by light entering the slit off center (for a 


detailed discussion on the effects of light entering 
the slit off-center see Maciejewski & Binney 2001). 
For example, consider the case in which the peak 
of line emission is on the blue side (left on Figure 1) 
of the NUC slit. Then the measured velocity will 
be blueshifted with respect to the true value. The 
same would happen for POSl. However in POS2 
the peak of the line emission will be on the red 
side of the slit and the measured velocities will be 
shifted to the red with respect to the true value. 
This is not the case for our data. The possible ori¬ 
gin and implications of this blueshift are discussed 
in §5. 

3.2. Morphology 

The acquisition image, shown in Figure 1, has 
a field-of-view of 5" x 5". In Figure 6 we show the 
inner 30" x 30" and 5" x 5" of the WFPC2 and 
NICS images. Finally in Figure 7 we plot the color 
maps. 

The acquisition image (Figure 1) shows a com¬ 
pact but resolved (FWHM~ 072) bright central 
feature superimposed to the central region of the 
bulge. A second bright feature is present ~ 074 
SW of the nucleus. Figure 8 represents the radial 
light profile from the STIS acquisition image. As 
already clear from the acquisition image, the light 
profile indicates the presence of a bright central 
feature. This feature is spatially extended as can 
be seen from a comparison with the light profile 
of an unresolved source (NGC 4051) obtained with 
the same instrumental setup. This nuclear feature 
is likely to be a star cluster; a photometrically dis¬ 
tinct star cluster is often present in the dynamical 
center of spiral galaxies of all Hubble types (e.g. 
Carollo, Stiavelli & Mack 1998; Carollo et al. 2002; 
Bbker et al. 2002, and references therein). 

The same central source is visible in all the 
WFPG2 images and an analysis of its colors (Fig¬ 
ure 7) indicates that it is bluer than the sur¬ 
rounding regions but still redder than the galaxy 
disk and bulge. Its Vega magnitudes in the three 
WFPG2 filters are 19.10 (F450W), 18.05 (F606W) 
and 17.16 (F814W) with formal uncertainties of 
±0.02mag. These values, converted to Johnson 
magnitudes, become B=19.2, R=17.7, 1=17.1 
for the KOV and Sb spiral spectra and B=19.1, 
R=18.0, 1=17.15 for the AOV and Sc spiral spec¬ 
tra. The star cluster emission is very weak in the 
K band and not readily visible but its location is 




Fig. 6.— Overlay of the K band isophotes on the grayscales of the F450W (left) and F814W (right) WFPC2 
images. North is up and East is left. The bottom panels show an expanded view of the nuclear region. 


coincident with the location of the K band peak 
within the alignment uncertainties (~ 0'.'2, see 
Figure 6). 

The color images show that the inner few arc- 
seconds are redder than the galaxy disk and bulge 
and this could either be due to the presence of 
obscuration by dust or to a change in stellar pop¬ 
ulation. 


From several galaxy catalogues it can be in¬ 
ferred that the inclination of the large scale disk 
is ~ 20°. This is supported by the nearly circu¬ 
lar isophotes observed in the NICS and WFPC2 
images at large scales > 20" — 30". However, in 
the inner regions, a bar-like structure is visible in 
the K band image oriented ~ E-W. It is extended 
for ~ 15" i.e. ~ 1.4 kpc with isophotes symmetri- 
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Jin i^rr»r] 

Fig. 7.— 20" X 20" view of the color images of the central region of NGC 4041. The center of the images 
coincides with the location of the nucleus. North is up and East is left. From left to right: B — I, B — I 
degraded to NIGS resolution, B — K with B (from WFPC2/F450W) degraded to NIGS resolution. The dark 
regions have redder colors. 



log Radius farcsec] 



log Radius [arcsec] 


Fig. 8.— Radial light profile from the STIS acqui¬ 
sition image. The empty diamonds represent the 
radial light profile from a galaxy, NGG 4051, hav¬ 
ing a strong central unresolved source; the profile 
was rescaled to match the NGC 4041 one in the 
central pixel. 

cally twisted on both ends. The twisting is most 
likely due to spiral arms just outside (and appar¬ 
ently connecting to) the bar. These give the bar 


Fig. 9.— Azimuthally averaged light profiles for 
the WFPC2 and NIGS images. All images have 
been rescaled to match the external light profile 
of the F814W image. The scaling factors applied 
are indicated in the figure and directly show the 
color in the outer regions where the profiles match. 
“NIGS res.” indicates that the F814W image has 
been degraded to the NIGS resolution. 
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the appearance of being larger than it is. It is 
well known (see, e.g. Shaw et al. 1995; Knapen et 
al. 1995) that strong concentrations of young stars 
can dominate the K-band. This likely happens in 
the region where spiral arms emanate from the 
bar in NGC 4041: blue-star complexes are seen in 
color images there, particularly at the west end 
(Figure 6). The bar itself is rather short with 
semi-major axis length at most ~ 3”, and a po¬ 
sition angle of about 60°. It is a weak bar, given 
the fatness of the observed structure and that it 
is lacking the pair of dust lanes characteristic for 
strong bars (Athanassoula 1992). Instead, a mini¬ 
spiral structure is seen in the STIS ACQ image 
(Figure 1), extending at least out to 3" from the 
nucleus, i.e. throughout the entire IR bar. Such 
a morphology indicates that we shouldn’t expect 
strong departures from the circular motion in the 
gas flow. In fact, the CO velocity field observed by 
Sakamoto et al. (1999, see their Figure If) shows 
no such deviations down to the 4" beam size of 
the observations. The PA of the kinematical line 
of nodes is ~ 85°, close to the PA of the bar. 

3.2.1. Light Profiles 

The canonical way to derive light profiles from 
images is to fit ellipses to the isophotes and to 
measure encircled fluxes. However, inspection of 
Figures 1 and 6 indicates the presence of small 
scale structures invalidating the ellipse fits of the 
WFPC2 and STIS images. The NICS image is 
much smoother but the ellipticity of its inner 
isophotes is caused by the presence of a bar and 
spiral arms and not by geometrical projection ef¬ 
fects. Since the inclination of the large scale disk is 
only ~ 20°, one expects almost circular isophotes, 
and it is reasonable to obtain the light profiles 
from azimuthal averages instead of canonical el¬ 
lipse fitting. 

In Figure 9 we plot the azimuthally averaged 
light profiles extracted from the WFPC2 and 
NICS images. The light profiles have all been 
rescaled to match the external one of the F814W 
image. The profiles all match well beyond r > 10" 
and the colors of that region can be estimated as 
B-I ^ 1.7, i? - / ~ 0.4 and / - a: ~ 2.1. These 
colors are consistent with other spiral galaxies: 
for instance de Jong (1996) has shown that the 
spirals of this Hubble type on average have inte¬ 
grated colors B — I = 1.8 ± 0.2 (observed 1.7) and 


H-a: = 3.5 ±0.3 (3.8). 

The red feature observed in the color maps 
around the nucleus results in a flattening of the 
light profiles which increases at decreasing wave¬ 
length. From the figure one can immediately 
estimate the color excess or reddening for each 
pair of bands with respect to the external points. 
Roughly E{B - /) ~ 0.8, E{R - /) ~ 0.4 and 
E{I — K) ~ 0.5. Assuming the galactic extinc¬ 
tion law by Cardelli, Clayton, & Mathis (1989) 
one finds that with Ay ~ 1.2 mag one can approx¬ 
imately match the observed color excesses (Ay ~ 
1.2 implies E{B - 7) ~ 1.03, E{R - J) ~ 0.32 and 
E{I — K) ~ 0.44). Thus the red feature observed 
in the color maps is likely to be due to extinction 
by dust with average value of Ay ~ 1.2 mag. 

4. Model fitting 

4.1. The stellar luminosity density 

In order to estimate the stellar contribution to 
the gravitational potential in the nuclear region, 
we have to derive the stellar luminosity density 
from the observed surface brightness distribution. 
This inversion is not unique if the gravitational po¬ 
tential does not have a spherical symmetry. Here, 
we assume that the gravitational potential is an 
oblate spheroid (e.g., van der Marel & van den 
Bosch 1998). We further assume that the prin¬ 
cipal plane of the potential has the same inclina¬ 
tion as the large scale disk {i ~ 20, LED A and 
de Vaucouleurs et al. 1991). The knowledge of 
f, combined with the observed axial ratio of the 
isophotes, should directly provide the axial ra¬ 
tio of an axisymmetric light distribution, but the 
presence of a bar-like structure invalidates this ap¬ 
proach since the flattening of the isophotes cannot 
be ascribed to simple projection effects. Therefore 
we consider two extreme cases: a spherical and a 
disk light distribution. The model light profiles 
are computed taking into account the finite spa¬ 
tial resolution and pixel size of the detector. A 
detailed description of the relevant formulas for 
the inversion and fitting procedure is presented in 
Appendix A. 

The light profiles derived from HST images 
have better spatial resolution than the one de¬ 
rived from ground-based NICS image. However, as 
shown above, the optical HST images suffer from 
Ay ~ 1.2 mag extinction in the nuclear region. It 
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Table 1 

Best fit parameters of the stellar light density. 
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pV 

b 
' b 

(3* 

^ a 

Pb 

mb ’’ 

a 

P 

X?ed 


I band (WFPC2/F814W) 

1.0 

1023 

0.27 

4.2 

0.73 

12.6 

0.9 

1.6 

10.8 

0.1 

951 

0.31 

5.0 

5.8 

12.6 

0.9 

1.6 

10.5 

K band (NICS) 

1.0 

910 

0.27'’ 

4.2^= 

0.44 

28.0 

1.4 

3.0 

9.8 

0.1 

609 

0.3P 

5.0^= 

1.60 

45.9 

1.4 

6.6 

7.7 

’’In units of M© pc 

assuming T = 

: 1. 




’’In units of arcsec. 

'’Kept fixed at the value from the previous fit. 


is also well known that K band light is a better 
tracer of mass. Therefore it is useful to use both 
the HST and ground based light profiles in order 
to infer the mass profiles. The star cluster is prob¬ 
ably washed out by the lower spatial resolution in 
the NICS image, thus, in order to account for its 
presence one can determine its geometrical param¬ 
eters from the fit of the HST light profiles and use 
them in the fitting the NICS image. 

In principle, the WFPC2 images could be red¬ 
dening corrected by using the color maps in Figure 
7 to derive the E{B — V) map. In order to do this, 
one might assume that the colors in the nuclear re¬ 
gion are constant and equal to those beyond 10" 
(see Figure 9). For a more detailed description of 
this procedure see for instance Marconi et al. 2000. 
However the cluster would be forced to have the 
same intrinsic color as the galactic disk, which is 
probably not realistic. Moreover its shape would 
be affected by PSF differences among the images. 
Finally it is not clear if the color differences are 
due entirely to reddening. We therefore decided 
not to apply any reddening correction and the con¬ 
sequences of this will be discussed in §5. 

We first consider the light profile derived from 
WFPC2/F814W because it is the least sensitive 
to reddening among the HST images. We then fit 


a model light distribution composed of the central 
star cluster and the more extended component. 
We assume that the star cluster luminosity density 
is spherically symmetric with a density law: 

( 1 ) 

The extended component has a functional form of 
the type 

p(m)=pj^\ ( 2 ) 

\mbJ y \mbJ J 

where jq^ corresponds to the 

radius in the spherical case and q is the axial ratio 
of the mass distribution. This functional form of 
the extended component is a reasonable descrip¬ 
tion of the central part of the galaxy where the 
bulge dominates. It is sufficient for our purposes 
because we are interested in modeling only the in¬ 
ner r < 5". 

In Equation 2, g = 1 for a spherical light dis¬ 
tribution, and g = 0.1 for a disk-like one. These 
are the two extreme cases. The best fit param¬ 
eters for the HST light profile are shown in the 
upper part of Table 1. The left panel in Figure 10 
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Fig. 10.— Left: fit of the light profile obtained from the WFPC2/F814W data. Right: fit of the NICS/K 
light profile with the central star cluster and the extended component (the geometrical parameters of the 
star cluster were determined in the previous fit). In both cases the lower panel shows the fit residuals. The 
dashed line in the right panel shows the K band fit obtained by fixing the cluster normalization to a value 
which is 0.6 dex larger than the best fit value. The empty squares show the corresponding residuals. 


shows the fit in the spherical case which is visually 
indistinguishable from the disk case. 

Using the geometrical parameters of the clus¬ 
ter determined from the fit to the HST data (r^ 
and (3*) we fit the NICS profile, both the spherical 
and disk case. The results are shown in Table 1. 
The right panels in Figure 10 show the fit in the 
spherical case which, as before, is visually indis¬ 
tinguishable from the disk-like case. 

The values of the reduced reported in Ta¬ 
ble 1 are much larger than the expected value of 
1 and the main reason can be found in the sys¬ 
tematic deviation of the residual around 3" (Fig¬ 
ure 10) where the data are systematically lower 
than the model. This systematic deviation, whose 
maximum value is ~ 0.03 dex, i.e. 7%, is present 
in both the optical and near-IR light profile and is 
caused by the presence of the bar-like feature. 

The fits to the light profiles allow an estimate 
of the I and K magnitudes of the star cluster. 
With the assumed distance, the cluster luminosity 
is lO"^-^ ^Qwi where wl refers to the F814W band 
(in this band the Sun has = 4.16). In the K 

band, the cluster has lO”^’^ spherical 

and 10’^-^ ^qk cases ( Mq^ = 3.286). 


These results suggest that the cluster could be 
blue, with estimated color J — RT ~ 0.9 (or 0.5 
in the disk case). The estimate of the K mag of 
the cluster is however uncertain due to the fact 
that the cluster is unresolved in the ground based 
observations. In order to set a limit to the amount 
of the cluster contribution to the K band light pro¬ 
file, in Figure 10 we also plot the results of the fit 
when the cluster normalization has been fixed to 
a value 0.6 dex larger (i.e. the cluster is 1.5 mag¬ 
nitudes brighter) than the best fit value (dashed 
line). The fit is worse, as clearly indicated by the 
residuals, and we consider this a firm upper limit 
for the cluster contribution to the K band light 
profile. This implies that the observed color of 
the cluster is I — RT < 2.4. If we assume that the 
cluster is subjected to Ay = 1.2 mag of extinction 
(§3.2.1), the upper limit changes to I — K < 2.0. 

Figure 11 shows the circular velocities in the 
principal plane of the potential expected from the 
density distributions determined from the fit of 
the NICS and HST data. In order to match the 
different rotation curves the circular velocities are 
plotted for T/f = 0.5 or T/ = 2.3, the values de¬ 
rived in §4.2 from fitting the rotation curves. It 
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Fig. 11.— Circular velocities in the principal plane 
of the potential for the density distributions deter¬ 
mined from the fit of the NICS and HST data. The 
circular velocities are derived assuming a spheri¬ 
cally symmetric mass distribution and are plotted 
for Tk = 0.5 and T/ = 2.3. The dotted line rep¬ 
resent the extreme case in which the cluster nor¬ 
malization in the K band has been fixed to -1-0.6 
dex of the best fit value (see text). 

is clear that the rotation curves are very similar 
beyond r > 0.5". However they differ at smaller 
radii because of the different contribution of the 
star cluster to the total mass budget. This differ¬ 
ence is due to the fact that the mass-to-light ratio 
has been assumed constant in each band. 

4.2. Kinematics 

In order to model the gas kinematical data (ve¬ 
locities and widths measured along the slit) we se¬ 
lect the simplest possible approach and we assume 
that the ionized gas is circularly rotating in a thin 
disk located in the principal plane of the galaxy 
potential. The latter assumption is not needed 
if the galaxy potential has a spherical symmetry. 
We assume that the disk is not pressure supported 
and we neglect all hydrodynamical effects. Thus 
the disk motion is completely determined by the 
gravitational potential which is made of two com¬ 
ponents: one is stellar and is completely deter¬ 



Offset from NUC center [arcsec] 

Fig. 12.— Contours of the modelled emission line 
flux distribution of the nuclear gas disk with over¬ 
plotted the slit positions (vertical dashed lines), 
the position of the kinematic center (filled square), 
and the line of nodes derived from fitting the kine¬ 
matic data (solid straight line). The dotted line 
is the line of nodes of the extended material. The 
ellipse indicates the nuclear gas disk outer limit 
for an inclination of 20° with respect to the line of 
sight. 

mined by the mass distribution, derived in the 
previous section, and by its mass-to-light ratio 
T. The other one comes from a dark mass con¬ 
centration (the black hole), spatially unresolved 
at HST-I-STIS resolution and defined by its total 
mass Mbh- This is the same approach which has 
been followed in all previous gas kinematical stud¬ 
ies of circumnuclear gas disks (e.g. Macchetto et 
al. 1997; van der Marel & van den Bosch 1998; 
Barth et al. 2001; Marconi et al. 2001). 

In principle, in order to constrain the BH mass, 
one should compare the emission line profile pre¬ 
dicted by the model with the observed one. Even 
in the simple potential described above, the line 
profiles can be very complicated, with multiple 
peaks (Maciejewski & Binney 2001). Neverthe¬ 
less, as shown in Figure 4, the line profiles are 
well represented by single gaussians (after sub¬ 
tracting the “blue” component) and therefore we 
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Fig. 13.— The assumed flux distribution com¬ 
pared with the data after folding with the tele¬ 
scope and instrument responses. The model values 
are connected by straight lines in order to guide 
the eye. 

compare the moments of the line profiles: the av¬ 
erage velocity (v) and velocity dispersion cr de¬ 
fined as cr^ = {v'^) — (u)^. In order to be com¬ 
pared with the observations, the model {v) and 
a are computed taking into account the spatial 
resolution of HST/STIS and the size of the aper¬ 
tures over which the observed quantities are mea¬ 
sured. The formulae used to compute velocities, 
widths and line fluxes and the details of their 
derivation, numerical computation and model fit¬ 
ting are described in Appendix B. To derive the 
observed (v) and cr, our approach is to parame¬ 
terize the observed spectra by fitting them with 
as many gaussian components as required by the 
data. Then, components which are not obviously 
coming from circularly rotating gas can be dis¬ 
carded (e.g. Winge et al. 1998) and the average 
velocities and line widths can be computed from 
the remaining ones (this is trivial if, as in the 


present case, only a single gaussian component 
is left). With this approach, the observed and 
model parameters, are compared. The advantages 
of our approach are twofold. Firstly, much com¬ 
putational time is saved in computing just (v) and 
a instead of the whole line profile. Secondly, mod¬ 
eling the observed emission line profiles in detail 
requires very high S/N data and it is seldom pos¬ 
sible to obtain such high quality data with HST 
except for a few galaxies (e.g. M87, Macchetto et 
al. 1997). Even with the high S/N available from 
8m class telescopes the matching of the line pro¬ 
files remains a significant problem (e.g. Cen A, 
Marconi et al. 2001). 

Thus, the model {v) and cr depend on E, the in¬ 
trinsic surface brightness distribution of the emis¬ 
sion lines, and on the following parameters: 

• So, the position of the kinematical center 
along the slit with respect to the position 
of the continuum peak; 

• 6, the distance between the reference slit cen¬ 
ter of the NUC slit and the kinematical cen¬ 
ter; 

• i, the inclination of the rotating disk; 

• 6, the angle between the disk line of nodes 
and the slits; 

• Vijys, the systemic velocity of the disk; 

• T, the mass-to-light ratio of the stellar pop¬ 
ulation; 

• Mbb., the BH mass. 

Not all of these parameters can be independently 
determined by fitting the observations. 

It can be inferred from the equations in Ap¬ 
pendix B that Mbh, T and i are directly coupled. 
This fact has already been discussed in detail by 
Macchetto et al. (1997) but here we briefly indi¬ 
cate how this coupling can be used to pose a lower 
limit on the disk inclination i. Considering the 
case where no BH is present, the amplitude of the 
rotation curves depends linearly on T*^'®sinz and 
its value can be determined by fitting the obser¬ 
vations. Therefore T will increase with decreasing 
i. An upper limit on T can be certainly set by 
the largest reasonable value for a very old stellar 
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Fig. 14.— Best fit standard model of the observed rotation curves (solid line) compared with the data. The 
dotted line is the best fit model without a black hole. The model values are connected by straight lines in 
order to guide the eye. Note that points from external and nuclear regions are not connected because they 
are kinematically decoupled. 

The right panel is a zoom on the nuclear disk region. The plotted model uses the mass density distribution 
derived from the K band light profile with the assumption of spherical symmetry. 


population (e.g. Bruzual & Chariot 1993; Maras- 
ton 2000) and this will immediately fix the lower 
limit for i. If a BH is present, T°-®sinz can still 
be determined from the kinematical data obtained 
beyond the BH sphere of influence. Of course, the 
underlying assumption is that both T and i are 
constant. 

4-2.1. The intrinsic surface brightness distribu¬ 
tion of emission lines 

The intrinsic surface brightness distribution of 
emission lines S is an important ingredient in the 
model computations because it is the weight in the 
averaging of the observed quantities. The ideal 
situation would be to have an emission line im¬ 
age with spatial resolution higher than the spectra. 


Not having that, Barth et al. (2001) have used a 
deconvolved HST image to model their HST data. 
This helped them to match the microstructure of 
the velocity field at large radii, but it is important 
to remember that the information about the BH 
mass comes predominantly from within the central 
unresolved source of unknown luminosity distribu¬ 
tion. Our approach is to attempt to match the 
emission line fluxes as observed in the spectra (or 
in emission line images) together with a variety 
of synthetic realizations of this unknown central 
emissivity distribution. We use this to estimate 
its impact on the derived value of Mbh- Thus, 
we assume simple analytical forms of the line flux 
distribution which well match the observed one 
along the slit, after folding with the telescope and 
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instrument responses as described above. 

It is clear from the right panel in Figure 3 that 
the flux distribution observed at the three slit po¬ 
sitions is very complex and cannot be described by 
a single radially symmetric component in the disk 
plane. Therefore we choose the approach of mod¬ 
eling directly the flux distribution on the plane 
of the sky by exploring various simple analytical 
forms of the line flux distribution. We take the in¬ 
trinsic light distribution at point x, y in the plane 
of the sky in the NUC slit reference frame {x is the 
position across the NUC slit while y is the position 
along it) to be defined as 

I{x,y)=Y,Io^f^( — ) (3) 

\Toi J 

where fi is a circularly symmetric function of char¬ 
acteristic radius Voi and weighting factor loi ■ Each 
component function fi is centered at (xoi,yoi), 
and Xi is the radial distance from this location 
Xi = [(x - Xoi)'^ + {y - VoiYf '^- This is an ap¬ 
proach similar to that of Barth et al. (2001) ex¬ 
cept that we use a synthetic realization of the line 
flux map. Although each bright patch is assumed 
to be circularly symmetric, its material is in cir¬ 
cular rotation about the galactic nucleus. Hence 
patches are constantly shearing rather than mov¬ 
ing as coherent units. 

The free parameters characterizing the flux dis¬ 
tribution are therefore (xoi, j/oi), loi and Xoi and 
are chosen in order to match the observed flux 
distribution along the slit after folding the model 
with telescope and instrument. For example, 
one of the functional forms used is the combina¬ 
tion of 4 exponentials and a constant baseline i.e. 
fiixi/xoi) = exp(ri/roi). Figure 12 shows the con¬ 
tours of this line flux distribution which matches 
the observed profile along the slit after convolving 
with the instrumental response (see Figure 13). 

4-2.2. The standaxd appxoach 

The standard approach followed so far in gas 
kinematical analysis is to assume that (i) gas disks 
around black holes are not warped i.e. they have 
the same line of nodes and inclinations as the more 
extended components, and (ii) the stellar popula¬ 
tion has a constant mass-to-light ratio with radius 
(e.g. van der Marel & van den Bosch 1998; Barth 
et al. 2001). In the present case, the blueshift of 


the inner disk (§3.1.2) indicates that the standard 
approach must generalized to allow for the kine¬ 
matical decoupling between inner and large scale 
disks. 

We use the emission line flux distribution de¬ 
rived in the previous section and we fix the incli¬ 
nation to i=20°, i.e. the inclination of the galactic 
disk. The free parameters of the fit are then Mbh, 
T, So, b, 9, 14ys and At4ys, the velocity shift of 
the extended component with respect to the nu¬ 
clear one. We perform the fit using the mass den¬ 
sity profiles derived for the I and K band, both 
in the spherical and disk cases. The results of 
the fit are shown in Table 2. Statistical errors at 
the 2a level on log Mbh and logT are ±0.2 and 
±0.1, respectively, thus the derived BH mass is 
Mbh = lto'7 X Mq. Figure 14 shows the best 
fit model (solid line) obtained with the mass den¬ 
sity distribution derived from the K band light 
profile with the assumption of spherical symme¬ 
try. The dotted line is the best fit model without 
a black hole. The left panel of Figure 16 shows the 
fit of the NUC data from the model with the mass 
distribution derived from the I band. Figures 14 
and 16, indicate that a model with no BH cannot 
account for the observed nuclear rotation curve, 
producing a velocity gradient which is shallower 
than observed. Note that the position angle of 
the line of nodes, a free parameter of the fit, is the 
same as the one inferred from the large scale CO 
velocity map by Sakamoto et al. (1999). This sup¬ 
ports the assumption that the nuclear disk is the 
continuation at small scales of the galactic disk. 
However the fit confirms that the nuclear disk is 
blueshifted by ~ 10 km s“^ with respect to the 
extended component and this argues against the 
assumption in previous sentence. The fit is greatly 
improved if a velocity shift of 8 km s“^ is allowed 
for the POS2 data. This velocity shift is the con¬ 
sequence of an error on the absolute wavelength 
calibration of the POS2 data and is well within 
the expected STIS performances (Leitherer et al. 
2001). 

The value of T in the K band derived from 
the fit varies between 0.2 (disk light distribu¬ 
tion) and 0.5 (spherical light distribution). This 
range of values is in good agreement with the typ¬ 
ical K-band mass-to-light ratios of spiral bulges 
(Moriondo, Giovanardi & Hunt 1998). Similarly 
the value of T in the I band ranges between 2.2 and 
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Table 2 

Model fit results in the standard approach. 


Band 

(a) 

q 

(b) 

log Mbh 
( c) 

T 

(d) 

So 

(e) 

b 

(e) 

e 

(f) 

^vs 

(g) 

ARsys 

(g) 

Xred 

(h) 

K 

1.0 

7.04 

0.52 

-0.04 

0.02 

-45 

1212 

11 

2.32 

K 

1.0 

0.0* 

1.00 

-0.09 

0.06 

-58 

1211 

11 

2.88 

K 

0.1 

7.09 

0.22 

-0.07 

0.08 

-32 

1213 

11 

2.13 

I 

1.0 

6.86 

2.28 

-0.01 

0.00 

-42 

1210 

12 

1.86 

I 

1.0 

0.0* 

2.65 

-0.03 

0.03 

-45 

1210 

12 

2.09 

I 

0.1 

7.14 

1.65 

-0.03 

0.01 

-47 

1212 

10 

2.22 


^ Band from which the mass density profile was derived. 

^ Assumed axial ratio of the mass distribution. 

Log of BH mass in units of Mq. 

Mass-to-light ratio in used band. 

® Arcsec 
^ Degrees 
® km s“^ 

'' Reduced x^- 

* The parameter was fixed to this value. 
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Table 3 

Model fit results in the alternative approach. 


Band 

q 

log Mbh 

T 

So 

b 

e 

i 

14ys 

Xred 

(a) 

(b) 

(c) 

(d) 

(e) 

(e) 

(f) 

(f) 

(g) 

(h) 


Extended Component 

K 

1.0 

0.0* 

0.48 

0.0* 

0.0* 

-43* 

20* 

1224 

2.7 

K 

0.1 

0.0* 

0.31 

0.0* 

0.0* 

-43* 


1224 

2.7 

I 

1.0 

0.0* 

2.29 

0.0* 

0.0* 

-43* 

20* 

1224 

2.5 

I 

0.1 

0.0* 

1.38 

0.0* 

0.0* 

-43* 

20* 

1224 

2.5 

Nuclear Disk 

K 

1.0 

<6.8 

1.29 

0.0 

-0.05 

-40 

20* 

1204 

0.66 

K 

1.0 

<6.6 

0.48^ 

^ -0.01 

-0.04 

-39 

35 

1205 

0.67 

K 

0.1 

<6.8 

0.88 

-0.01 

-0.05 

-40 

20* 

1205 

0.69 

I 

1.0 

<6.8 

4.05 

-0.01 

-0.04 

-41 

20* 

1206 

0.60 

I 

0.1 

<6.8 

3.77 

0.00 

-0.04 

-41 

20* 

1205 

0.58 


^ Band from which the mass density profile was derived. 

^ Assumed axial ratio of the mass distribution. 

Log of BH mass in units of Mq. 

Mass-to-light ratio in used band. 

® Arcsec 
^ Degrees 
s km s“^ 

'' Reduced x^- 

* The parameter was fixed to this value. 
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Fig. 15.— Best fit alternative model of the observed rotation curves compared with the data. The solid line 
is the best fit model obtained by fixing i and varying T while the dotted line represent the case in which i 
has been varied and T kept fixed. The model values are connected by straight lines in order to guide the 
eye. Note that points from external and nuclear regions are not connected because they are kinematically 
decoupled. The right panel is a zoom on the nuclear disk region. The plotted model uses the mass density 
distribution derived from the K band light profile with the assumption of spherical symmetry. 


3.6 in agreement with measurements within the in¬ 
ner 2 kpc of spiral galaxies (Giovanelli & Haynes 
2002 ). 

4-2.3. Alternative model: the decoupled nuclear 
disk 

To date, everyone who has determined a central 
BH mass from gas kinematics assumed that the 
disks are unwarped, i.e. coplanar at all radii, and 
that the stellar mass-to-light ratio is constant with 
radius. 

The high surface brightness and velocity offset 
of the nuclear disk with respect to the extended 
component presented in §3 indicate that, at least 
for NGC 4041, this assumption could be untrue. 
The nuclear and large scale disks might be char¬ 


acterized by different geometrical properties like 
position angle of the line of nodes, inclination 
and systemic velocity. Recently Cappellari et al. 
(2002) have shown a discrepancy between the BH 
mass in IC 1459 determined from gas kinematical 
(Verdoes Klein et al. 2000) and stellar dynamical 
models. The authors propose that a possible solu¬ 
tion to this discrepancy could be an error on the 
assumed gas disk inclination. 

Therefore, in the alternative approach pre¬ 
sented here, we first fit the extended component 
data in order to determine the mass-to-light ratio 
T of the stellar component to be used in the fit of 
the nuclear data. For the extended component we 
assume a constant line flux distribution and deter¬ 
mine T in the two extreme cases of the spherical 
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Fig. 16.— Best fit models of the observed rotation curves computed using the mass densities derived from 
the I band light profiles. The model values are connected by straight lines in order to guide the eye. The left 
panel refers to the standard model (the dotted line is the model without a BH) while the right panel refers 
to the alternative model. 


and disk-like light distribution (i.e. with an axial 
ratio of g = 0.1). 

From the CO velocity map by Sakamoto et al. 
(1999) the PA of the kinematic line-of-nodes is 
PA= 86°, i.e. almost exactly East-West. Since 
the PA of the slit in the STIS observations is 43°, 
this means that the angle between the slit and 
the disk line of nodes is 9 = —43° (for the defini¬ 
tion of 9 see Appendix B). We also assume that 
the inclination of the disk is z=20°, as specified in 
the previous section. The CO velocity map also 
clearly indicates that the large scale velocity field 
is circularly symmetric. 

The fitted circular rotation curves of the ex¬ 
tended gas are shown in Figure 15. The derived 
T for the spherical or disk geometry are shown in 
Table 3 for the mass density profiles derived from 
both the K and I band light. In our derivation 


we allowed for a constant velocity shift between 
the data points in the off-nuclear slit positions 
and the nuclear one. This is caused by the uncer¬ 
tainty in the absolute wavelength calibration, and 
is very small, -1 km s“^ and 8 km s“^ for POSl 
and POS2 respectively. 

Similarly to the results of the standard model, 
mass-to-light ratios in the K and I band are con¬ 
sistent with typical values for the central regions 
of spiral galaxies. 

We now fit the nuclear data (r < 0'.'4) allowing 
for values of T or z different from those of the 
extended component. Results of the fit are shown 
in Figure 15 and in Table 3. The right panel of 
Figure 16 shows the fit of the NUC data from the 
model with the mass distribution derived from the 
I band. 

In all these cases the BH mass has been set to 
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Fig. 17.— Statistical upper limits on the BH mass 
in the alternative approach for the cases in which 
the mass density has been derived by the I and K 
band light profiles. 



Position along the slit [arcsec] 

Fig. 18.— Observed and model FWHMs at the 
NUC position. The solid line is the case without a 
BH (alternative model) and the dashed line is with 
a 10^ M© BH (standard model). The assumed 
intrinsic velocity dispersion is 35 km s“^. 


zero and the goodness of the fit suggest that no 
BH is needed to fit the data with this approach. 
We estimate the upper limit on the BH mass by a 
1 parameter variation, i.e. we assign a fixed value 
of Mbh and we perform the fit. When we have 
>= 1,4,9 we have reached the 1,2,3 a upper 
limit. As an example, in Figure 17, we show the 
plot corresponding to the K and I band spherical 
cases. There the 3a upper limit is 10® ® Mq. 

From Table 3 it can be seen that the position of 
the nucleus along and across the slit is independent 
of the assumed mass distribution. The inclination 
varies, but like T, is just a scaling factor. Indeed 
varying T produces the same results, as shown in 
the same table. The angle between the slit and 
the line of nodes is very well determined in the 
range —39 — —41°. This means that the PA of 
the kinematical line of nodes is 43° (the PA of 
the slit) minus —40° i.e. ~ 83° consistent with the 
PA of the kinematical line of nodes on large scales. 
The systemic velocity is 1204—1206 km s“^ and is 
definitely blueshifted with respect to the extended 
rotation where it is 1224 km s“^. 

As pointed out by Barth et al. (2001) the 
FWHM of the lines might be a worrisome issue 
in the sense that FWHM much larger than in¬ 
strumental might indicate motions which could 
deviate from circular (apart from the broadening 
of the line due to unresolved rotation around a 
point mass). In Figure 18 we plot the observed 
FWHMs and compare them with the expected 
one in the framework of the alternative model 
assuming an intrinsic velocity dispersion of only 
35 km s“^. Also the dashed line represents the 
case with a 10^ Mq BH showing that, in this case, 
the FWHM cannot pose good constraints on the 
BH mass. 

5. Discussion 

5.1. How the BH mass determination is af¬ 
fected by assumptions. 

In §4.2 we have shown that the inclination 
of the nuclear disk, f, and the mass-to-light ra¬ 
tio of the stellar population, T, play a critical 
role in detecting the presence of a nuclear BH. 
The standard model, in which i and T are the 
same for the nuclear disk and the more extended 
galactic disk, requires the presence of a BH with 
Mbh = lla? ^ (§4-2.2). Conversely, al- 
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Fig. 19.— From top to bottom: I — K,\ogM/L k 
and log M/Lj vs time for a single stellar popula¬ 
tion (instantaneous burst) with solar metallicity 
(models by Bruzual & Chariot 1993, updated in 
1996). The solid and dashed line are for a Salpeter 
(1955) and a Scalo (1986) IMF respectively. In the 
upper panel, the solid horizontal line represents 
the observed upper limit on the I — K color. The 
dotted line is the same upper limit after correction 
for Ay = 1.2 mag. In the central and lower panel, 
the horizontal dotted lines indicates the range of 
mass-to-light ratios allowed by the observations. 


lowing either T or i to vary, the kinematical data 
can be fit without a BH (the so-called alternative 
model, §4.2.3) and Mbh < 6 x 10® Mq. These 
results show a caveat in gas kinematical searches, 
where it is always assumed that i and T are con¬ 
stant at all radii. 

To distinguish between the standard and alter¬ 
native models it is necessary to find some extra 
physical constraints. In particular the value of 
mass-to-light ratio T determined from the obser¬ 
vations should be consistent with realistic stellar 


populations (see §4.2). Little can be said about i 
which in contrast with the other geometrical pa¬ 
rameter of the disk, 9, is poorly constrained by 
the rotation curves. Usually one can infer an es¬ 
timate of i from the morphology of the disk ob¬ 
served in emission line images, especially in E and 
SO galaxies (e.g. van der Marel & van den Bosch 
1998; Barth et al. 2001). The root of the prob¬ 
lem lies in the coupling between Mbh, "f and i 
that arises because they are all essentially scaling 
factors on the amplitude of the rotation curves. 

In the standard approach (Table 2) Mbh does 
not depend on the assumed mass density profile 
(either from the K or I band). The angle 6, which 
is a free parameter of the fit, is quite well con¬ 
strained around —40°, i.e. the PA of the kinemat¬ 
ical line of nodes is the same as the one at large 
scales derived from the 2D CO map by Sakamoto 
et al. (1999). The constancy of 6, which is not an 
intrinsic parameter of the system but only depends 
on the galaxy orientation in the plane of the sky, 
suggests that i {Ai between the nuclear and galac¬ 
tic disk is an intrinsic parameter of the system) 
should also be constant. No similar argument can 
be placed on T which can then be allowed to vary 
and it is then possible to fit the kinematical data 
without the presence of a BH. This requires that 
T in the nuclear region is a factor of ~ 2 larger 
than that in the extended disk. 

To verify if this T is consistent with the ob¬ 
served colors, we plot in Figure 19 I — K, Tk and 
T/ as a function of time from the burst for a single 
stellar population experiencing an instantaneous 
burst of star formation. We have used the models 
by Bruzual & Chariot (1993), updated in 1996, we 
considered solar metallicities, Salpeter (1955) and 
Scalo (1986) IMF (solid and dashed lines in the fig¬ 
ure), and used the theoretical stellar libraries. The 
solid line in the upper panel marks the upper limit 
we placed on the I — K color of the stellar clus¬ 
ter (§4.1); the dotted line marks the same upper 
limit but dereddened for Ay = 1.2 mag (§3.2.1). 
Finally, the dotted lines in the central and lower 
panels shows the T required to fit the data in the 
alternative model both for the spherical and disk¬ 
like distribution of stars. The Ts implied by the 
analysis are thus consistent with a very old nu¬ 
clear star cluster, with age ~ 10 Gyr. This is ap¬ 
parently in contrast with the high Ha equivalent 
width measured in the nuclear spectrum, ~ 70 A, 
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which seems to require the presence of very young 
stars. A possibility is that the Ha emission is not 
due to young stars but to a low luminosity AGN 
and, indeed, the total Ha luminosity in the nu¬ 
clear region is L(Ha) = 4.9 x lO"* L© correspond¬ 
ing to the low end of the A (Ha) distribution for 
AGNs (Ho, Filippenko, & Sargent 1997). To prove 
that the emission lines are excited by an AGN, 
one could use the so-called diagnostic line ratios 
[NH]/Ha and [OHIJ/H/3 (e.g. Osterbrock 1989). 
Ground based spectra with 1" resolution (Axon 
et al. in preparation) indicate that this object is 
on the boundary between HII regions and LIN- 
ERs with [NH]/Ha ~ 1 and [OHI]/H/3 < 0.5. 
Unfortunately the spectra coverage of our HST 
data is limited to 6500-7000 A region but the the 
[Nil]/Ha ratio at O'.'l resolution remains similar 
to that observed from the ground suggesting that 
the same could be true for the [OHI]/H/3. 

If on the other hand the star cluster is young, 
it can fully account for the Ha emission. As¬ 
suming case B recombination, the A (Ha) lumi¬ 
nosity implies an ionizing photon rate of (5(H) ~ 
2 X I0®° s“^ (e.g. Osterbrock 1989), which corre¬ 
sponds to ~ 2 — 20 O stars (Panagia 1973). When 
confronted with stellar population synthesis model 
(Leitherer et al. 1999) the observed Ha equivalent 
width indicates an age of ~ 10^ yr. In this case, 
as shown in Figure 19, the I band mass-to-light 
ratio is more than a factor ten lower and the ht 
of the kinematical data requires a BH even in the 
alternative model. However when using the stel¬ 
lar mass density profile derived from the K band, 
the data can still be explained by stars alone. The 
reason is that the contribution of the star cluster 
to the total light density prohle is small, as can be 
seen from the rotation curves in Figure 11. Note 
that in order to explain the high [Nil]/Ha ratio 
with photoionization from stars one needs to allow 
for higher-than-solar N abundances in the nucleus 
as has been suggested many times in the past since 
the work by Burbidge & Burbidge (1962). 

In order to proceed further and firmly es¬ 
tablish the age of the star cluster, one needs 
HST/NICMOS infrared images of the nuclear re¬ 
gion. The optical to near-IR colors of the star 
cluster constrain its age and mass-to-light ratio. 
Moreover with HST/STIS blue spectra one can 
measure the [OHIJ/H/3 ratio, in order to estimate 
the AGN contribution, and verify the existence of 


deep H Balmer lines in absorption expected if the 
cluster is young. 

The other assumptions behind the modeling are 
the shape of the stellar density profile (spherical 
vs disk-like, I vs K band, effects of reddening), the 
intrinsic flux distribution of the emission lines and, 
finally, the assumption of circular motions. 

The stellar density profile depends on q, the 
intrinsic axis ratio of the extended mass distri¬ 
bution, and on the relative importance between 
the nuclear star cluster and the extended stellar 
component. As shown in the previous section, the 
four cases examined there cover the most extreme 
cases: spherical or disk-like distribution of the ex¬ 
tended component {q = 1 and g = 0.1 respec¬ 
tively), stellar cluster which gives a small contri¬ 
bution to the total stellar density distribution (K 
band) or which dominates over the nuclear disk 
size i.e. r < 0"4 (I band). From Tables 2 and 3 
we can conclude that the value of the BH mass 
or its upper limit do not depend on the assumed 
stellar density profile. From the same table it can 
be seen that the position of the nucleus along and 
across the slit is also independent of the assumed 
mass distribution. In §4.1 we decided not to ap¬ 
ply a reddening correction to the light profiles be¬ 
cause of the uncertainties involved. This redden¬ 
ing correction would have affected the values of 
the mass-to-light ratios derived from the rotation 
curves (see §4.2) whose accurate measurement is 
not the aim of this paper. Roughly, assuming an 
average Ay = 1.2 mag (§3.2.1), the I light would 
increase by ~ 70% and T/ would decrease by the 
same amount (-0.2 in log). The K light would 
similarly increase only by ~ 10% with a similar 
decrease in (-0.04 in log). These corrections 
would have a marginal effect for the conclusions on 
the cluster age derived in Figure 19. Finally the 
reddening correction effects on the shape of the 
stellar light/mass density profiles, would be negli¬ 
gible for the estimate of the BH mass upper limit. 
Indeed, the same value is obtained regardless of 
the use of mass density prohles derived from the I 
or K band light which are very differently affected 
by reddening. 

Another important issue is the influence of the 
intrinsic line flux distribution on Mbh- We have 
tried several other different realizations of the ob¬ 
served flux distribution by varying both the num¬ 
ber of components and the functional form of 
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single components (exponential, gaussian, power 
law). We have found that accurately describing 
the observed flux distribution is important. Mod¬ 
els that did not do this produced fits significantly 
worse than those described in the previous section. 
When the fits of fluxes and velocities are accept¬ 
able, Mbh or its upper limit do not depend on the 
assumed flux distribution. 

The assumption of circular motions is proba¬ 
bly one of the most critical ones. Circular mo¬ 
tions in gaseous disks are expected due to the dis¬ 
sipative nature of the gas, however anisotropies in 
the stellar potential rapidly lead to non-circular 
streaming (e.g., Athanassoula 1992). The effect of 
non-circular motions would be that of “distorting” 
the rotation curves and increasing the intrinsic line 
width of the gas. Therefore, obtaining very good 
fits (Xred < 1) of the rotation curves in the al¬ 
ternative model strengthens our confidence on the 
assumption of circular motions. This assumption 
is also supported by the small intrinsic line width 
(ct ~ 35 km s“^) required by the observations. In¬ 
deed cr/wcirc can be used to quantify the effects 
of turbulence or non-circular motions on the BH 
mass determination (e.g. Verdoes Klein et al. 2000; 
Barth et al. 2001). For the best fitting models 
f/ucirc is always less than 0.4 (provided that there 
is s BH with a mass close to the upper limit in the 
case of the alternative models). This is similar to 
what Barth et al. (2001) found in the case of their 
best fit model of NGC 3245. There the effect of 
the asymmetric drift correction, a possible way to 
include non-circular motions in the analysis (e.g. 
Binney & Tremaine 1987), is that to increase the 
estimate of the BH mass by just 10%. Apart from 
these qualitative arguments, the presence of sig¬ 
nificant non-circular motions cannot be excluded 
and it is not possible at the moment to quantify 
their effects on the method we have adopted and 
described here. 

Non-circular motions are certainly present in 
the nuclear region of NGC 4041 as indicated by the 
presence of the blue wing, but these have been sin¬ 
gled out in the deblending procedure as in Winge 
et al. (1998). 

5.2. The blueshift of the nuclear disk 

A proposed explanation for the blueshift ob¬ 
served in the nuclear disk is that the star cluster 
is oscillating across the galactic plane. In this pic¬ 


ture, the cluster is bound to the galaxy, is old and 
has a large T ensuring that it is massive enough 
and not subjected to tidal disruption. 

Any velocity component perpendicular to the 
disk will give a large contribution to the observed 
velocity because the galaxy is almost face-on. If 
the cluster is oscillating across the galactic plane, 
the observed blueshift may be translated into a 
velocity modulus of similar magnitude. Thus, the 
cluster velocity (~ 10 — 20 km s“^) is smaller than 
the rotational velocity (see Figure 11) and the 
cluster is bound to the galaxy. If the cluster is 
very massive, as in the alternative model, then the 
gravitational potential over the star cluster size is 
dominated by the star cluster itself and it is not 
subjected to tidal disruption. Also in this picture 
the gaseous disk is completely dominated by the 
gravitational potential of the star cluster. The to¬ 
tal cluster mass within r < 073 as derived from 
the I-band mass density profile (the one in which 
the cluster dominates) is ~ 6 x 10^ M© and this 
gives an average density of ~ 200 M©/ pc“^, not 
an unusually high value, since it can be found in 
galactic globular clusters (Pryor & Meylan 2002). 

In the picture in which the cluster is young 
(from the fit of the K band we have I — K = 0.9, 
i.e. ~ 10^ yr from Figure 19), the stability against 
tidal disruption could pose a problem, since the 
mass in the nuclear region is completely domi¬ 
nated by the bulge stars. However, this problem 
could be solved if the lifetime of the cluster against 
tidal disruption is less than or equal to the cluster 
age. 

5.3. Mbh - galaxy correlations 

It is now clear that a large fraction of hot 
spheroids contain a BH and, moreover, it seems 
that the hole mass is proportional to the mass 
(or luminosity) of the host spheroid. Quantita¬ 
tively, Men/Msph « 0.001 (e.g. Merritt & Fer- 
rarese 2001). This relation is still controversial, 
however, both because the sensitivity of published 
searches is correlated with bulge luminosity, and 
because there is substantial scatter in Mbh at 
fixed Mgph. Recently Ferrarese & Merritt (2000) 
and Gebhardt et al. (2000) have shown that a 
tighter correlation holds between the BH mass 
and the velocity dispersion of the bulge. The two 
groups however find two different slopes of the cor¬ 
relation, Mbh cx ct* and Mbh cx respectively 
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(Tremaine et al. 2002). Clearly, any correlation of 
black hole and spheroid properties would have im¬ 
portant implications for theories of galaxy forma¬ 
tion in general, and bulge formation in particular. 
Indeed, several authors have shown that the slope 
of the Mbh-ct* correlation yields information on 
the formation process. If the formation process is 
self-regulated, i.e. if the BH growth is limited by 
radiation pressure, Mbh oc ct* (Silk & Rees 1999). 
Conversely if the growth is regulated by ambient 
conditions Mbh oc af (Adams, Graff & Richstone 
2000; Cavaliere & Vittorini 2001). The uncertain¬ 
ties on the slope of the Mbh-ct* correlation do not 
allow one to distinguish between the two cases. To 
solve this problem more BH mass measurements 
are needed in the low mass range (I0®-10^ Mq) 
where spiral galaxies are expected to fit. 

These correlations are also very important to 
estimate BH masses quickly and easily instead 
from very complex dynamical and kinematic mea¬ 
surements. Therefore the quantities involved in 
the correlations must be measured as carefully as 
possible in order to reduce the scatter to its intrin¬ 
sic value. For instance it has been shown that with 
careful estimates of the bulge luminosity the Mbh- 
Lgph correlation has the same scatter as the Mbh- 
CT* correlation, in contrast with previous claims 
(McLure & Dunlop 2002). 

Given the low value of the BH mass in NGC 
4041, Mbh < 2 x 10^ Mq, it is worthwhile to 
verify the relation of this galaxy with the pro¬ 
posed correlations. The B magnitude from the 
RC3 catalogue is 11.9 becoming 11.8 after extinc¬ 
tion correction (see NED). The morphological type 
is Sbc/Sc and T = 4.0 ± 0.3. From Simien & 
de Vaucouleurs (1985) the bulge to total luminos¬ 
ity ratio is ~ 0.16 resulting in Am = 2. Thus 
the bulge magnitude is 13.8. The adopted dis¬ 
tance, 19.5 Mpc, corresponds to a distance mod¬ 
ulus of 31.5 thus the absolute bulge magnitude is 
— 17.7 corresponding to 1.8x10® Lhq. The best 
fit of the Mbh — Ae.sph correlation gives Mbh = 
0 . 8 xl 08 (LB,sph/ 10 iOLBe)io 8 (Kormendy & Geb- 
hardt 2001). Thus the expected BH mass in NGG 
4041 would be 1.2x10^ Mq, in agreement with 
the BH mass estimate or upper limit. The best 
fit of the Mbh — cr* correlation gives Mbh = 
1.3X 10®(cr*/200 km s“^)^ ® (Tremaine et al. 2002) 
or Mbh = 1-4 x 10®(cr*/200 km s“^)^ ® (Merritt 
& Ferrarese 2001). Using INTEGRAL/WYFFOS 


at the WHT, we have recently measured the stel¬ 
lar velocity dispersion in the central 2" of NGG 
4041 (Batcheldor et ah, in preparation). With 
CT* = 95 ± 5 km s“^, the expected BH masses are 
then 7x10® M© and 4x10® M©, both consistent 
with the result from this paper. 

Since the main goal of our project is to de¬ 
termine whether or not spirals do in fact follow 
the MBH-Asph and Mbh-ct* relations, it is impor¬ 
tant to observe also objects which, a-priori, are ex¬ 
pected to have marginally or non detectable BHs. 
Indeed, even if the present measurement is only 
an upper limit, this is still useful in ruling out the 
presence of unusually massive central BHs in late 
type spiral galaxies. 

6. Conclusions 

We presented HST/STIS spectra of the Sbc spi¬ 
ral galaxy NGG 4041 which were used to map the 
velocity field of the gas in its nuclear region. We 
detected the presence of a compact (r ~ 074 ~ 
40 pc), high surface brightness, circularly rotating 
nuclear disk cospatial with a nuclear star cluster. 
This disk is characterized by a rotation curve with 
a peak to peak amplitude of ~ 40 km s“^ and is 
systematically blueshifted by ~ 10 — 20 km s“^ 
with respect to the galaxy systemic velocity. 

We have analyzed the kinematical data assum¬ 
ing that the stellar mass-to-light ratio is con¬ 
stant with radius and that the gaseous disk is not 
warped, having the same inclination as the large 
scale galactic disk. We have found that, in order 
to reproduce the observed rotation curve, a dark 
point mass of (IIq’t) ^ 10^ ^0 is needed, very 
likely a supermassive BH. 

However the observed blueshift suggests the 
possibility that the nuclear disk could be dynami¬ 
cally independent. Following this line of reasoning 
we have relaxed the standard assumptions varying 
the stellar mass-to-light and the disk inclination. 
We have found that the kinematical data can be 
accounted for by the stellar mass provided that ei¬ 
ther the mass-to-light ratio is increased by a factor 
of ~ 2 or the inclination is allowed to vary. This 
model resulted in a 3cr upper limit of 6 x 10® M© 
on the mass of any nuclear black hole. 

Gombining the results from the standard and 
alternative models, the present data only allow us 
to set an upper limit of 2 x 10^ M© to the mass of 
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the nuclear BH. 

If this upper limit is taken in conjunction with 
an estimated bulge B magnitude of —17.7 and 
with a central stellar velocity dispersion of ~ 
95 km s“^, the putative black hole in NGC 4041 
is not inconsistent with both the MeH-T^sph and 
the Mbh-ct* correlations. 
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A. Deriving the luminosity density of the stars from the observed surface brightness profile 


The stellar luminosity density can be inferred by inverting the observed observed surface brightness 
profiles. Following van der Marel & van den Bosch (1998) we assume an oblate spheroidal density distribution 
which we parameterize as: 

m is defined as m? = x'^ + y'^ + j, where xyz is a reference system with the xy plane corresponding 
to the principal plane of the potential, q is the intrinsic axial ratio of the mass distribution. If T is the 
mass-to-light ratio, the observed surface brightness distribution is given by^"^: 

1 /•+“ 


where the integration is performed along the line of sight, ft can be shown that 


1 q 1"^°° p{m^)d'm? 

47rT q' J ^,2 ^ {jn? — 


(A3) 


where = x'^ + y'^ jq'^ and x'y' is a reference system on the sky, with the x' axis aligned along the 
apparent major axis, q' is the observed axial ratio of the isophotes which is related to the intrinsic axial 
ratio of the mass distribution by q'^ = cos^ i + q^ sin^ i, where i is the inclination of the line of sight (z = 90° 
is the edge-on case). The observed surface brightness results from the convolution of S with the point spread 
function P of the system (i.e. telescope and optics) and the detector pixels: 


/'X-\-A.X pY-\-XY ■ ^+oo ^+oo 

^app(A, r)=/ / / / T^tmeix, y)P{x'- x,y'- y)dxdy 

Jx-AX JY-AY U-ao J -oo 


dx'dy' 

IKxKy 


(A4) 


where X, Y is the centre of an aperture with size 2AA x 2AF. The integration on dx'dy' can be directly 
carried out on the PSF P which is described by a sum of gaussians: 


P{x,y) 



(A5) 


with Nk = ^nafki. Then Eapp is given by 

f+°° r+°° dxdv 

Sapp(A,y) = y J Y{x,y)truelC{X-AX-x,X + AX-x-,Y-AY-y,Y + AY-y) {A6) 

where /C is the convolution kernel: 

1 w , 

/C(Ai ,X2;Y,,Y2) = —Y. 2^^? i 

i—1 

where 8 is the complementary error function. The integration is then carried out numerically using the 
Gauss Legendre approximation. 
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Note the l/47r factor. This derives from the fact that p/V is a density (luminosity per unit volume) and E is a surface brightness 
(luminosity per unit area per unit solid angle). 
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The best fitting model is determined by minimizing the reduced Xied defined as: 


Xred 


J_ / log Sj - log Sm (r,, (5r,; Pi, ■. ■, Pm) \ ^ 
S log Si J 


(A8) 

(A9) 


where i = 1,..., N indicates a data point with surface brightness Si ± JSi at radius ri. 

Sri;pi,... ,Pm) is the model surface brightness, averaged over radii ri — 6ri < r < n + Svi, which is 
a function of m free parameters pi, ■ ■ ■ ,Pm- = N — m is the number of degrees of freedom. The x^ed i® 
minimized to determine the m free parameters using the downhill simplex algorithm by Press et al. (1992). 

When the mass density has been determined as described above, the circular velocity in the principal 
plane is given by the relation (e.g. Binney & Tremaine 1987): 


V;^{r) = [ 

Jo 


p(m?)m?dm 


(AlO) 


where e is the eccentricity related to g by = 1 — e^. 
defined by m < rrio is 

/ rtio 


Similarly the mass enclosed within the homoeoid 
p{rri^)m^dm (All) 


B. The rotation curve model 

The velocity field along the line of sight, v, can be easily computed in the case of a circularly rotating thin 
disk. Let XY be a reference frame in the plane of the sky with Y axis fixed along the direction given by the 
slit position (hereafter called the slit reference frame). Consider a reference frame X^Yd in the sky with the 
Xd axis aligned along the disk line of nodes and the origin coincident with the disk centre (hereafter called 
the “disk reference frame”, see Figure 20). The origin of XY is chosen in such a way that the disk center 
has coordinates x = b and t/ = 0 in the slit reference frame. Then a given disk point P with coordinates 
{x,y) in the slit reference frame has coordinates {xd,yd) in the disk reference frame given by 


Xd = {x — b)sin9 + ycosO 
yd = —{x — b)cos9 + ysm0 


(Bl) 


If the disk has an inclination angle i (i = 0 in the face-on case), then P is at the disk radius r given by 


r^ = x‘'d+(—) 

V cos z/ 


The circular velocity of P, in the case of a spherical mass distribution, is then given by 


V'c(r) = 


dd) 

dr 


2 ^GM(r)^2 


(B2) 


(B3) 


where M{r) is the enclosed mass at radius r, a constant value Mbh in case of a point mass. In the case of 
an oblate spheroidal mass distribution, Vc(r) is given by Equation AlO. 

The velocity component along the line of sight is finally given by 

V = 14ys - Vc(r) sinz^ = V^ys - (GMbh)°'^ (B4) 

with the latter expression describing the simple case of a point mass Mbh- 
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Fig. 20.— Geometry of the disk. 


This velocity field v has to be convolved with the instrumental response in order to simulate the observed 
quantities. Below we show how can this be done for any velocity field v. 

Consider again the reference frame XY in the plane of the sky defined above. The light distribution can 
be written as 

^{x',y',v) = I{x',y')(j){v-v), (B5) 

where I{x',y') is the total intensity at the point {x',y'), and (l){v — v) is the intrinsic line profile centered 
at the velocity along the line of sight v{x',y'). After passing through the telescope, the light gets convolved 
with the instrumental PSF P{x — x',y — y'), and at point {x,y) in the focal plane the light distribution is 
described by 


[x,y,v) = 


dx'dy' I(x', y')4‘{'v — v)P{x — x',y — y'). 


As already defined, the slit in the focal plane is aligned along the Y axis and let its center cross the X axis 
at xq. Let the slit width be 2Aa;. Past the slit, the light falls into the detector plane, where the spatial 
coordinate x and the velocity v are combined into one single detector coordinate w = v + k{x — xo), which we 
identify with the observed velocity (we denote it w here in contrast to the intrinsic velocity denoted by v). 
Here, the coefficient k is given by iiAw/Ay where /r, the anamorphic magnification, accounts for the different 
scales on the dispersion and slit directions. In the case of STIS, the scale along the slit is 0'.'05071 and along 
dispersion is 0705477, thus y = 1.080 (Bowers & Baum 1998). The light distribution in the detector plane, 
d>{w,y), is calculated by ousting n, and integrating the light contribution across the slit 


rXo-\-Ax 

d>{w,y) = / dx ^{x,y,w — k{x — Xo)) 

J Xq — ^X 
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nXo + Ax p p + oo 

/ dx'dy'(j)[w - v{x',y') - k{x - xo)]I{x',y') P{x - x',y - y'). 

Jxq—Ax J J—oo 


(B7) 


Note the component k{x — Xg) in the velocity profile, which is the “spurious” velocity for light entering off 
from the slit center. Properties of thus defined two-dimensional light distribution with the spurious velocity 
shift were investigated in detail by Maciejewski & Binney (2001). 

The detector integrates over finite pixel sizes therefore we calculate the expected line fluxes, average 
velocities and widths for the line profile that is obtained by integrating the light distribution on the detector 
plane d^(w,y) over the width 2Ay of the pixel along the slit, and convolving it with the shape of the 
pixel in the dispersion direction: a top-hat of width 2Aw. Thus the expected line profile in the detector is 


= 


pVj+Ay 


pw-\-Aw 


lyj-Ay 


dy 


> w — Aw 


dw d!{w', y), 


(B8) 


where the pixel has the coordinate yj along the slit. Note that the observed intensities are measured at 
discrete values of w corresponding to the pixel centers. In order to calculate the expected line fluxes, average 
velocities and widths, one has to evaluate moments of 


(•-1-00 


dw 


^j{w) = 


rVj+Ay 


fXo+Ax 


/ dy dx 

yj—Ay Jxq — Ax 


+ 00 /*+oo 

tA.J -n / 


dx'dy' V 


pw-\-Aw 

w" / dw'(j)(w' — Wg) (B9) 

J w — Aw 


where we used abbreviations V = I(x',y')P(x — x',y — y') and wg = v{x',y') + k(x — xg). The last two 
integrals can be simplified by inverting the order of integration: 


/ -poo pw-\-LAw p f-roo 

dw ic" / dw'(j){w' — Wg) = / / dw dw' w'''(j){w' — Wg)H{w — w') = 

-OC Jw — Aw J J—OQ 

(w + Aw)"~^^ — (w — Aw)"~^^ 


c+oo 


dw (j){w — Wg) 


(^w+Aw 


/ w — Aw 


C + oo 


dw 4>{w — Wg)- 


n+\ 


where H{w — w') is 1 for \w — w'\ < Aw and 0 otherwise. This leads to the following formulae for the 
moments of 


c + oo 


'^j{w)dw = 2 Aw 


fXo+Ax 


fVj+Ay 


I dx dy 

xq — Ax Jyj—Ay 


+ 00 


dx'dy'V 


/ + 00 _ pXg + Ax 

wj{w)dw = 2Aw / dx 

-oo Jxq — Ax ^ 

(*+oo __ nXQ+Ax nyj-\-Ay 

w"^ ^j{w)dw = 2Aw / dx dy 

— OO Jxq — Ax Jyj—Ay 


rvj+^y 


dy 


Vj-^V 

*+oo 


+ 00 


dx'dy'{i 


dx'dy'Wg V 


{Aw)'^ 


'■)V 


(BIO) 

(Bll) 

(B12) 


To obtain them, we assumed that the intrinsic velocity profile 4'{v) is bound and symmetric, i.e. 


/ -l-oo p-roo /*-hOO 

(j>{v)dv = 1 / v(l>{v)dv = 0 / v‘^(l){v)dv = 

-OO J —OO J —oo 

where cr^ is the intrinsic velocity dispersion. 
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Equations BIO, Bll and B12 can be used to compute the expected line fluxes, average velocities and 
widths. For example, in the most simple case of a constant line intensity / = const, and velocity held 
V = const, one can write: 

_ _ 

(B13) 

f-oo 3 3 

Then the expected velocity dispersion, given by (t| = (v'j) — (vj)'^, is, understandably, larger than the 
intrinsic velocity dispersion cr. This is due to the convolution with the pixel size and slit width, which add 
quadratically. However, note that while Aw and a enter the integral (B12) as constants, wq is a linear 
combination of the spurious velocity shift, and the intrinsic velocity. These two contributions can cancel, 
resulting in the expected line profile being broadened by the pixel size only, and not by the width of the slit. 
This implies that the wide slit can probe the disk on the scale of the pixel size rather than the slit width. 
Maciejewski & Binney (2001) explore consequences of this finding. 

In order to compute the model given by Equations BIO, Bll, and B12, we create a grid in x and y with 
sampling given by a-psp/n. Here, ctpsf is the spatial r.m.s. of the point spread function (PSF) and n is 
the subsampling factor. We have verified that the optimal subsampling factor used is n = 3 since larger 
values do not produce any appreciable differences in the final results. The PSF used is the one generated 
by TinyTim (V6.0, Krist & Hook 1999) at 6700A. Convolution with the PSF is done using the Fast Fourier 
Transform algorithm Press et al. (1992). Following Barth et al. (2001) we have also introduced the CCD 
scattering function Leitherer et al. (2001) but, as already noticed by them, it does not have any appreciable 
effect in the final results. 

We compare models to the observed spectrum, which is essentially an array of intensities d'y, after the 
observed line profile is derived in the following way: to the sequence of intensities 'I'y for a given row j along 
the slit, we fit a baseline, and a continuous analytical function which we interpret as the observed 

equivalent of the expected line profile dtj(w) (eq.BS). The best fitting model is determined by minimizing 
the reduced Xred defined as: 



-Ass 


fc=i i=i 


^kj - {Vj)k{pi, ■ ■ ■ ,Pm) 
Svkj 


IFfcj (yVj)kiPl^ • • • j Pm') 


where the index fc = 1,3 indicates the slit position, and j = l,Nk counts pixels along the slit. Here, the 
characteristics of the model are as follows. The velocity in the row along the slit {vj)k, and the FWHM 
of the velocity profile {Wj)k, are calculated directly from equations (B13), now clearly for variable intensity 
and velocity held. They both are functions of m free parameters pi, ■ ■ ■ ,Pm, which are determined by Xred 
minimization. The FWHM is calculated from the expected velocity dispersion aj after assuming a Gaussian 
line profile. The observed velocities {vki ± Svki), and velocity dispersions {Wki ± SWki) are also derived from 
equations (B13), but after ^j(rc) has been replaced by ^j(rc)°’^® defined above. ‘^^k — m is the 

number of degrees of freedom. 

The Xred minimized to determine the m free parameters using the downhill simplex algorithm by Press 
et al. (1992). In order to apply statistical methods when Xred i® much larger than 1, we follow Barth et al. 
(2001) and rescale errors as: 

i6v'ki)^=dvl+6V^ (B15) 

where 6V is a “systematic” error determined such that the resulting Xc of tho “best” model is 1. The fits 
presented in this paper have Xred ~ f ff*® error rescaling was not performed. 
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